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Abstract 

We measure the pressure and energy density of two flavor QCD in a wide 

range of quark masses and temperatures. The pressure is obtained from an 

integral over the average plaquette or (ipip). We measure the QCD (3 function, 

including the anomalous dimension of the quark mass, in new Monte Carlo 

simulations and from results in the literature. We use it to find the interaction 

measure, e — 3p, yielding non-perturbative values for both the energy density 

e and the pressure p. 
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I. INTRODUCTION 



We expect that at high temperatures strong interactions will enter a new phase, the 
quark-gluon plasma (QGP), which is believed to have existed in the extremely high temper- 
atures microseconds after the big bang. Heavy-ion collision experiments at the Brookhaven 
AGS and CERN are currently trying to recreate the QGP. In order to prove its existence in 
the aftermath of a heavy-ion collision and to understand the dynamics of the QGP in the 
early universe one needs as input, among other things, the equation of state for the system. 
Since the phase transition occurs in a regime of strong gauge coupling, a non-perturbative 
method is called for. Lattice calculations provide such a method. However, at currently 
practical lattice spacings the operator formalism usually used in such calculations for the 
equation of state is not necessarily non-perturbative, since it requires the knowledge of the 
asymmetry coefficients, or Karsch coefficients [I]. These are currently known only perturba- 
tively [2]. These asymmetry coefficients are short distance quantities defined at the scale of 
the lattice spacing a, so that if a could be made small enough (and temporal size Nf large 
enough to keep the temperature fixed) the perturbative coefficients could be used accurately 
even though the temperature would remain at a typical QCD length scale. However, the 
use of perturbative values for the asymmetry coefficients leads to distortions in the equation 
of state at the bare couplings used today. An effort can be made to non-perturbatively 
measure the asymmetry coefficients [3]. In practice that has turned out to be difficult [4]. 

The integral method does not require the knowledge of these coefficients. It was first 
used in the context of lattice QCD to calculate the interface tension by S. Huang et al. 
[5] and later modified for the bulk pressure of pure gauge QCD by J. Engels et al. [6]. 
The disadvantage of the integral method is that for the pressure at a single temperature and 
quark mass, a number of different simulations are required in order to provide the integrand. 

We have done an extensive survey of the gauge coupling and quark mass plane for two 
fiavor QCD, allowing us to measure the non-perturbative pressure by integration. Using 
data from the literature for the p and tt meson masses, we calculate the non-perturbative fi 



2 



function to compute the interaction measure and hence the energy density. 

Ahhough these simulations avoid one of the problems of using QCD simulations on small 
lattices, namely the use of perturbation theory, the problems of scaling violation, or non- 
constant ratios of physical lengths, of flavor symmetry breaking with the Kogut-Susskind 
quarks, and of effects on the thermodynamics of replacing integrals over the momenta by 
sums over discrete Matsubara frequencies remain. 

In Sec. II we present the formalism for calculation of thermodynamic quantities and the 
/3-function. Section III details our simulations and the results for the energy and pressure. 

II. THEORY 

A Euclidean X Nf lattice with periodic boundary conditions has a temperature T and 
volume V given by 

V = Ny, 

l/T = N,a , (I) 

where a is the lattice spacing. The form of the partition function Z for QCD with rij 
fermions with Kogut-Susskind (KS) discretization is 

Z = I [dU^n,,)]exp{ {6/g')S, + {rij/A) Trlog[am, +P] } , (2) 

where the gauge fields Ui^n,^) are SU(3) matrices at the link (n, fi) (at site n, to the direction 
= 0, ...,3), is the gauge coupling and aniq the bare quark mass in lattice units. The 

gauge action 

•S'fl = ^Re J2 ^rUa{n,fi,v) , (3) 

is a function of f/n(n, ;/), the path ordered product of link matrices around the elementary 
plaquette at site n in the fip plane. The covariant derivative contains the Kogut-Susskind 
phases. For nj = 2, the simulation is performed using the standard refreshed molecular 
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dynamics algorithm [7] . This involves integration of an equation of motion with a nonzero 
step, and a resulting error in physical averages. As it turns out, this step size error must be 
handled with care. 



A. Thermodynamics 

Thermodynamic variables are derivatives of the partition function Z defined in Eq. (2). 
In particular, the pressure p and energy density e are given by 

eV=-^^ (4) 

d{\IT) ^ ' 

and 

V_ _ dlogZ 
T~ dV • 

For large, homogeneous systems the free energy is proportional to the volume: 

9k)gZ 

~dv~ ■ 

Thus, the free energy density / can be connected to the pressure by 

pv , 7 -fv 

— = log Z = -jr- ■ (7) 

The free energy cannot be obtained from a single simulation, but its derivatives can be. In 
particular, 

— 1 Slog Z 

" 2N^Ntd{&lg^) ' 

where (□) is the average plaquette normalized to 3 for a lattice of unit matrices and 

These are relatively easy to measure in a simulation. If a series of runs is performed with 
different aniq and 6/(/^ values the pressure can be obtained by numerically integrating Eqs. 
(8) and (9). With the plaquette we obtain 



^(6/(/2, am,) = 2N^Nt [(□(G/^/'^, am,)) - (□(G/^/'^, am,)) ] <6/(/'2) , (10) 

-/ ^cold sym 

where (□(G/*/^))^^^ is the average plaquette from a symmetric (cold) system. The subtrac- 
tion of the symmetric value removes the divergent zero temperature pressure. The lower 
limit for the integration should be in a region where the difference between the zero and 
nonzero temperature plaquette expectation values is negligible. This removes the unknown 
constant introduced by the integration. In a completely analogous way we get from 

^;-(6//,am,) = iV3iVW I g\rn'^a)) - I g\w;^a)) \d{m'^a) . (11) 

1 ./cold sym 

or 

;^(am,) = r\U4^{m'^a)) - • (12) 

1 J cold sym 

At high temperatures, pjT'^ should approach a constant. This means that the integrand 
in Eq. 10 must approach zero at large 6/(/^. Figure 1 shows the A^t = 4 plaquettes and cold 
lattice plaquettes as a function of 6/(/^ at am, = 0.1. As required, the curves join at large 
6/(/^ values. 

Although the average plaquette curves coalesce, the difference between the spatial and 
temporal plaquettes approaches a constant as shown in Fig. 2, and as required by the 
operator formula for the entropy. This means that at very high temperatures the spatial 
and temporal plaquettes are shifted by equal amounts but in opposite directions from the 
zero temperature plaquette. 

B. Energy density 

To obtain the equation of state for QCD we also need the energy density e. It can 
be obtained non-perturbatively using the /3-function and the interaction measure /. The 
interaction measure / is 

IV eV pV 1 dlogZ dlogZ 
T T T Td{l/T) dV 
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FIG. 1. The Nt = 4 plaquettes (octagons) and the symmetric plaquettes (squares) as a function 



of at arriq = 0.1. 
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FIG. 2. The difference of spatial and temporal plaquettes as a function of at anig = 0.1. 
The diamond corresponds to 16"^ X 4 lattice, the squares to 12"^ X 4 and the octagons to 8"^ X 4 



lattices. 
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= y^^a -^^^^i ^ = 

oat oiog a 

= 2N;N.^m-(n)„J + N;N.^\{i.^t)-{i''t) 1 , (13) 

oiog a ^ oiog a ^ ' ^ ' sym 

where and at are the spatial and temporal lattice spacings. In the last step we have 
subtracted the zero temperature value. We have measured (□) and {^V'V') on the lattice and 
calculated the /3-function from mass spectrum data in the literature. 

Knowledge of the non-perturbative pressure and interaction measure allows us to com- 
pute other bulk quantities. The energy and entropy s become 

e = I + 3p (14) 
sT = I + 4:p . (15) 

Given the need for the /3 function to find the interaction measure from Eq. 13, or for the 
asymmetry coefficients to find the energy directly, it would be nice to compute the energy 
density by an integration similar to that used for the pressure. To the extent that we can 
make a small change in temperature by changing A^t, this is possible. Let us start from 
Equation (4) 

dloK Z A log Z 

eV -- 



d{l/T) A{Nta) 

N'^ /-amq 

= 7Wr^WV iv;[-(V'V)^, + (V'V') ] + iV,[(V'V)^ -(V'V') ' (16) 

(iV/ — iVtjaJcold ^ ' ^ 'sym ^ ' \ 'sym « 

where and are two different temporal extents. Of course, an analogous formula in 
terms of the plaquette also exists. 

Equation (16) can be given in a form 



e 



where we have taken as the inverse temperature the average [N^ + Nt)a/2. 

Ideally, we would like to take N^a and Nfa as close as possible to each other to increase the 
accuracy of the approximation of the derivative by the finite difference. This approximation 
gives a curve that is smoother than the real energy density. In fact, it gives the average e/T'* 
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in the temperature range [1/ {Nla),l / [Nta)). For linear regimes of p/T"^ it is exact. This is 
true, for example, in high temperatures, where the Boltzmann law is expected to be valid. 

Approximating the pressure at 1/T = [N^ + Nt)a/2 with the average of 1/T' = N^a and 
1/T" = Nftt systems, we can get a formula for the entropy s as well: 

s [NlNt{Nl + Nt)/2f /•»™.. 



[(V'V')^ - {i^i^)jdmy (18) 



The entropy does not need a vacuum subtraction; the symmetric average cancels out. 

In the applications of the equation of state, the sound velocity of the thermal system is 
an important quantity. Acoustic perturbations travel in the system with a speed c^: 

1 de 
c1 dp 

One has to take the derivative keeping the physical quark mass fixed, i.e., on the line of 
constant physics. Unfortunately, we know best only the variations of the energy density 
and pressure along lines of constant bare parameters. In order to measure the correct sound 
velocity one has to use the fi function to map the changes in bare parameters to physical 
changes of temperature and quark mass. 
Define 



4 

e = ea 



p = pa'^ 



I = la*, (201 



as the dimensionless quantities measured on the lattice. Then, the sound velocity becomes 



1 _ dl-AI^^ 
dp — Ap 



3+ .ja ^ (2i: 



where the differentials have to be taken along a path which keeps mj^/nip constant. There- 
fore, it is not enough to have the energy density and pressure as the function of bare pa- 
rameters; one needs the /3 function as well. We will return to this in future work. 



8 



C. The /3-function 



To change the lattice spacing keeping physical ratios fixed, we must adjust both relevant 
couplings in the action, 6/(/^ and am^, so our /3-function has two components: 



am„ 



\ dlog a ' Slog a 



(22) 



where a is the lattice spacing. This can be measured from and nip in units of the lattice 

spacing as functions of and aniq. In general, a change in the bare quantities changes 

both the lattice spacing a and the physical quark mass. In practice, to keep physics constant 

while changing a, the partial derivatives in Eq. 22 are taken at constant m^/mp. Clearly, the 

nucleon mass or any other physical mass could be substituted for mp, and in the continuum 

limit should give the same answer. In parctice, is special since it is uniquely sensitive to 

the quark mass. (In principle, other quantities sensitive to aniq could be used, such as the 

nucleon-delta mass splitting.) 

We have used tt and p masses from simulations with two fiavors of Kogut-Susskind 

fermions reported in the literature (see Table I) to calculate the /3 function over a wide 

range of coupling and quark mass. A 2 X 2 matrix M was formed by calculating derivatives 

from three or four points in coupling-quark mass space: 

(91n(am^) (91n(am^) 
(9(6/g2) (9(amq) 



d In(am^) 




d In(amp) 





(91n(amp) (91n(amp) 
d{Q/g'^) (9(amq) 





-d{%lg^)- 


X 






d (aniq) 



(23) 



The partial derivatives were found by linear interpolations among the points. The /3 func- 
tion was then calculated by inverting Eq. 23 with the changes on the left hand side set equal: 
d In(am^) = d In(amp) = 6. This keeps mj^/nip constant. We have also calculated the /3 



function at = 5.35 and 



0.1 from new simulations [9]. In these simulations, we 



varied space and time couplings anisotropically to measure the asymmetry coefficients. The 
asymmetry coefficients give the change in couplings as the space and time lattice spacings 
are adjusted independently. However, their sums give the symmetric change in the couplings 
when all lattice spacings are varied together, or the usual /3 function. Values for the nonper- 
turbative /3 function are given in Table II, and a plot of the renormalization group (RG) fiow 
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FIG. 3. Renormalization group flow in the gauge coupling-quark mass plane for two flavors of 
Kogut-Susskind fermions. The base of each line is indicated by an octagon and is where the (3 
function is evaluated. The end is indicated by a one standard deviation error ellipse. The inset 
shows this error ellipse for the /3 function at = 5.65 and arrig = 0.0175, the point indicated by 
the arrow. The length of each line and its error ellipse corresponds to a scale change of c? ln(a) = 0.1, 
so this is actually 0.1 times the error on the (3 function. The two lines at = 5.35 correspond 
to the (3 function calculated using the VT p and PV p masses. (The larger line is for the VT p.) 
The differences are due to flavor symmetry breaking from the large lattice spacing. 
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is shown in Fig. 3. The values of 6/(/^ and aniq quoted in Table II are averages of the points 
used in calculating the fi function. In Table II we also give the value of the quark mass 
component of the fi function minus the contribution from the classical dependence on the 
lattice spacing. The RG flow depicted in Fig. 3 shows how to approach the continuum limit. 
From Table II we see that the gauge piece of the fi function is roughly half the perturbative 
value at couplings for the range of 6/(/^ used in these simulations. 

The errors shown in Fig. 3 were calculated by forming a singular co variance matrix for the 
two components of the fi function for each of the hadron masses anih used in the calculation: 



These matrices were added together to obtain a nonsingular covariance matrix. The covari- 
ance matrices for the first two entries in Table I were calculated from a jackknife estimate. 
The covariance matrices are diagonalized, and the allowed variance is then given by an 
ellipse whose semimajor and semiminor axes are along the eigenvectors of the covariance 
matrix and whose lengths are given by the square roots of the corresponding eigenvalues. 
The ellipses shown in Fig. 3 were calculated for one standard deviation. The error ellipses 
in Fig. 3 are close to straight lines, which indicates the two components of the fi function 
are highly correlated. Correlations between the tt and p masses from the same simulations 
were not included in our analysis except for the point at 6/(/^ = 5.35. 

There are two systematic errors in our calculation of the fi function. The first comes 
from the linear approximation of the matrix of derivatives. However, Fig. 3 shows that 
the fi function is smooth over a wide range in coupling- quark mass space. This indicates 
the first order approximation is good even for these large changes. The second and most 
apparent systematic error comes from scaling violations in the masses (see 3). The point at 
6/(/^ = 5.35 and aniq = 0.1 was calculated using both the VT p and the PV p masses, and 
each gives a different answer for the fi function. (See Ref. [8] for a discussion of the meson 
operators.) The difference is attributable to the slopes of the masses as functions of 6/(/^ 
and aniq. Since the VT p and PV p are not degenerate at this point but should be in the 
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continuum limit, this is expected. At this couphng the slopes differ roughly by a factor of 
two, which leads to the large discrepancy in the /3 function. At larger 6/(/^ and smaller aniq 
the problem is alleviated since the p masses become degenerate to good accuracy. 

Since we are interested in the /3 function at many points in coupling - quark mass space 
to evaluate the interaction measure, a better method is to fit the spectrum data as a function 
of and aniq. Such a fit can also be used to transform functions of the bare quantities 
to functions of physical parameters like the temperature. Thus, we determine m^/mp and 
niptt as functions of aniq and The inverse function then yields the /3 function. The 

fitting functions are given in Table III and shown in Fig. 4 where we have used nip = 0.770 
GeV to convert to inverse lattice spacing. These fitting forms are ad hoc fits to the masses 
in the relevant parameter region, and do not have the correct asymptotic behavior for large 
or large aniq. We take the functional form of mj^/nip from chiral perturbation theory 
with coefficients that are polynomials in We fit the mass ratio in this case because we 

could not obtain a fit with reasonable for the pion mass alone. These fitting functions 
are compared to the simulation results for am^ and anip in Table I. For the region of 6/g^ 
relevant to the Nf = 4: thermal crossover, it can be seen that this is a good interpolating 
function for these masses. The results (not shown) are also fit reasonably well by a 
quadratic function of 6/g^ and aniq. In the continuum limit, the functions in Figs. 4(a) 
and 4(b) should converge, i.e., = m^2- However, at the couplings used in current lattice 
simulations, a pronounced breaking of the continuum SU[2) X SU[2) chiral symmetry is 
evident. 

The /3 function is determined by computing numerical derivatives of the bare parameters 
with respect to lattice spacing along lines of constant m^/mp, i.e., from the inverse of the 
functions shown in Fig. 4. The confidence levels in Table III are low, which may indicate 
that the errors on the spectrum data are underestimated. In any case, the fits should be 
considered as giving smooth interpolated mass values in the parameter regions where data 
exist. Fvidence that they work is the agreement between the resulting /3 function and that 
from the direct calculation as shown in Fig. 5. We note that there seems to be a discrepancy 
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at 6/(/^ = 5.35 and aniq = 0.1. This is where the /3 function was calculated from new 
simulations. These simulations show that the p mass is a much flatter function of the bare 
parameters. As the quark mass is reduced below 0.1, it steepens appreciably. Since we have 
no spectrum data between aniq = 0.1 and aniq = 0.05, our fits can not resolve this behavior. 
In Table II we also give the /3 function from the fitted spectrum at a few selected points, 
including extrapolations to zero quark mass. 

The errors on the /3 function from the fitted spectrum were calculated as described above 
for the direct method. That is, we determined a singular covariance matrix for each mass 
used in the fit at each point where the /3 function was evaluated and proceeded as before. 




0.5 1.0 1.5 2.0 1 2 

(a) 1/a (GeV) (b) l/a (GeV) 

FIG. 4. Spectrum fit showing contours of constant bare coupling and quark mass, (a) Using 
the mass of the Goldstone pion of the exact U(l) X f^(l) lattice chiral symmetry, (b) Using a 
non-Goldstone pion mass. The fitting functions are given in Table III. The dashed line is the 
physical value of m-,^ jmp. The inverse function gives the renormalization group fiow. The octagons 
indicate the points where zero temperature spectrum calculations were done. The solid lines are 
concurs of constant bare parameters 6/^^ and aniq with the values indicated in the graphs. 
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FIG. 5. Comparison of QCD nonperturbative (3 function calculated directly from spectrum data 
and from the lit to the spectrum data. The octagons are the direct calculation, and the diamonds 
are from the fit. The locations of the octagons are averages of the points in parameter space used 
to calculate the (3 function. The diamonds show where thermodynamic data were obtained for the 
equation of state. 



14 



III. SIMULATIONS 




FIG. 6. The bare parameters of our runs. The plot symbols indicate the step size used, where 
the cross is dt = 0.01, the squares 0.02, the octagons 0.03, the diamonds 0.04 and the plusses 0.05. 
Some of the points have runs with several step sizes. The plusses connected by the solid line show 
estimates for the phase transition or crossover region for Nt = 4 [18,19]. 



We performed simulations with the parameter values displayed in Fig. 6. At each point 
we ran both hot [N^ = 4) and cold [N^ = Ng) lattices. For < 5.45 we used A^ = 8 or 
12. For 5.45 < < 5.69 we used A, = 12, and for = 5.77 we used A, = 16. On the 
cold lattices we performed about 800 trajectories plus 100 warmup trajectories, on the hot 
runs about 1600 plus 100 warmup trajectories 

As mentioned in Section 2, these results are subject to step size errors. We found that 
these are not the same for the hot and cold runs. In cold lattices the effect was much more 
pronounced. Thus it is not safe to subtract the results of cold and hot lattices without making 
sure that the step size errors are under control. This increases the workload considerably. 

The system, with fixed step size dt is self contained, that is it corresponds to some 
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statistical system and the integrations along constant aniq and 6/(/^ should give consistent 
results also at non-zero dt, provided our integrations are numerically accurate. So, this self 
consistency, although not strictly speaking physical, can be used to test the accuracy of our 
integrations. Indeed at fixed step size, the integrations agreed within error bars. 

In our case the proper handling of step size errors was especially important in the in- 
tegration over 6/(/^. This is because the difference of the hot and cold plaquettes was in 
many cases of the same order as the step size error. In Fig. 7a we show the plaquette 
as a function of step size squared, which is the leading error in the R-algorithm [7]. The 
difference in the step size errors between hot and cold runs also has implications for the op- 
erator measurements of pressure and energy density where a similar zero point subtraction 
is needed. 
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FIG. 7. (a) The variation of the plaquette as a function of step size squared at Qjg^ = 5.45 and 
mqa = 0.025. The diamonds come from S"* and the octagons from 8"^ X 4 lattices. The cold system 
has a much more pronounced effect, (b) The same for tptp. The cold system again has a larger 
effect, but it is much less significant for the calculation of the pressure, since the shift in tptp due 
to the step size error is a smaller fraction of the difference between the hot and cold values. The 
squares and bursts are runs with = 12. 
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For t/jt/j the step size error was not as large relative to the difference in values for the 
hot and cold condensate, as is seen in Fig. 7. Thus, the error in the measurement of the 
pressure from the aniq integrations is smaller. However, for the smallest values, = 5.35 
and aniq = 0.025, the extrapolation of t/jt/j to dt = had to be taken into account: with 
step sizes 0.02, 0.03 and 0.04 the cold lattice V'V' was 0.3021(20), 0.3328(25) and 0.3867(33) 
respectively. Linear extrapolation in step size squared gave t/jt/j = 0.2731(28), significantly 
different from the smallest step size result. In contrast, the hot t/jt/j was almost the same at 
dt = 0.02 and 0.03. This was the worst case in the mass integrations. 

For the integrations the step size errors are hardest to handle. The effects at 

aniq = 0.1 were within error bars. At m^a = 0.025 there were two major effects. First, the 
cold plaquettes were smaller with larger step sizes, increasing the apparent difference between 
the hot and cold systems. This means, that at larger step sizes the pressure apparently 
becomes too large. Second, the position of the phase transition was shifted for the hot 
system towards smaller 6/ as step size is reduced. Therefore, the integrated pressure near 
the transition became smaller at larger step sizes, partially cancelling the first effect after 
integrating to large 6/g^. In all, the step size error caused the pressure to look steeper than 
it really is. In Fig. 8 we display the pressure for two different step sizes. To get the physical 
value we make a linear (in step size squared) extrapolation of the pressure to zero step size 
also shown in the figure. The extrapolation is in agreement with the values of pressure from 
the mass integrations. In our final results, when a point could be reached by two integration 
paths, we combined the results of the two paths with weights proportional to the inverse 
squares of their statistical errors. Roughly speaking, the effect of this on the pressure versus 
6/g^ is to use the integrations over am^, with smaller statistical errors, as "anchors" for the 
integrations over 6/g^. 

Figures 9 and 10 collect our results for the pressure with = 4. In Fig. 9 we have 
combined the results of integrating over 6/g^ and aniq. 

Having the results as a function of am^, we may wish to extrapolate to even smaller, 
physical quark masses. For small quark masses aniq the t/jt/j should go as 
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ho{6/g^) + hi{6/g^) ■ arrig + 0{amgf,T > 
y^tP= €0(6/5^) + ci(6//) • anig + 0{amgY, T = (25) 

where hi and Ci are independent of mass, with ho = for T > Tq- Using Eq. (12) one gets 
for the pressure 

pa^{amg)= {pqO^) + / [{ho - Cq) + {hi - Ci){am'g)]d{am'g) 

= (poa*) + (ho - CoKanig - m") + [(amj^ _ (^")2] + 0(am3), (26) 

where the lower hmit m" is smaU enough for Eq. (25) to be vahd and {poa"^) is the value 
of the pressure at that point. Therefore, a simulation is needed only down to a mass value 
where the behavior of Eq. (25) is established. 

However, if the transition is of second order, near the critical temperature the scaling 
of t/jt/j is governed by the critical exponent 6 and these formulas must be altered. Using 
Eq. (26) to extrapolate the curves in Eig. 10 results in the bursts in Eig. 9. 

It is interesting to notice that at high temperatures, the behavior of the pressure becomes 
[ho ^0, m" ^ ) 

pa'^{amg) = poa'^{mg = 0) - V'V'coid'T^ga'* (27) 

and its mass derivative is determined by the zero temperature t/jt/j. In physical units 

p{mg) = p{0) - (^tptfj^^^^mg (28) 

Thus, even in very high energy scales, the zero temperature subtraction produces a nonzero 
derivative of the pressure with respect to quark mass at anig = 0, in contrast to free quark 
behavior where this slope is zero. 

The energy density from the interaction measure and /3 function using Eqs. 13 and 14 is 
displayed in Eigs. 11 and 12. The /3 function was obtained from the fits to m^/mp and nip 
in Table III. 

The finite size error in the pressure comes from the assumption (Eq. 6) that the partition 
function scales with the volume. However, since most of the simulations were done in a region 
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where the finite size efi'ects in expectation values are small, we expect only small corrections. 
Even the integration through the crossover region is mostly at a large mass, where finite size 
effects are not so important. The correlation length is determined by the quark mass aniq 
rather than by the lattice size Ng. The integration is over the crossover region at a 

small mass, where one might find an effect. But since the and aniq integrations agree, 
the effect cannot be too large. 

The largest systematic error comes from the fact that = A lattices contain only a few 
of the Matsubara frequences in the partition function. One can estimate this by comparing 
the finite free field theory results to continuum results. The Boltzmann law for massless 
quarks in the continuum gives 

^ = —(16 + 21) = 12.1725. (29) 
r4 30^ ^ V ; 

For free Kogut-Susskind quarks on a 12"^ X 4 lattice at anii = 0.025, the energy is 

^(iV, = 4,iV. = 40) = 1.46^16 
^^^P^(iV, = 4,iV. = 40) = 1.51^21 (30) 

giving e/T'* = 18.1. This is close to the value measured by us. 

Fig. 4 shows how to translate our results to physical units (the coordinate axis can be 
changed to temperature by dividing by Nf). The lines of constant aniq tend up as 
increases since the physical quark mass becomes larger as the lattice spacing a is reduced at 
constant aniq. Fig. 4(a) shows that the most direct approach to the physical value of quark 
mass (or m^/mp), is to reduce aniq. However, Fig. 4(b) reminds us that at small 6/(/^, even 
though the ratio of the Goldstone pion mass to the rho mass is near the physical value, the 
other pions are much heavier and the fiavor symmetry is badly broken. 

In Fig. 13 we show the pressure in physical units. 

Finally, we look at the different contributions to the state variables from the fermionic 
and gluonic sectors. Since we have used the integration method for the pressure and did not 
directly measure the energy, we only look at the interaction measure. In Fig. 14 we show the 
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contributions from both terms in Eq. 13 separately for the two hnes of constant m^a. Both 
pieces are normahzed to zero at T = by vacuum subtraction. As the transition region is 
traversed, the gauge piece shoots up while the fermion piece rises more slowly and more or 
less levels off. At high temperature the interaction measure is expected to go to zero. It is 
interesting to see how the two contributions approach this limit. In the high temperature 
phase, (V'V') approaches zero and the fermion contribution is given by the vacuum subtraction 
(a constant) times the quark mass component of the /3 function which goes to zero in the 
continuum (T — )• oo) and chiral limits. Fig. 5 shows why this contribution is not falling 
off for niqa = 0.1. The quark mass component of the /3 function is actually increasing in 
this region of while for m^a = 0.025 it is decreasing. The gauge contribution behaves 
in just the opposite way. The gauge component of the /3 function goes to a constant while 
the average of the plaquette in the hot phase approaches the zero temperature value (see 
Fig. 1). 

In this study we have developed and tested methods for determining the equation of state 
for high temperature QCD, and presented results with a large lattice spacing a = IjAT. 
Simulations must be done at many values of 6/(/^ and am^, and extrapolations made to the 
physical values. Zero temperature simulations must be done to allow the divergent parts 
of the energy and pressure to be subtracted, to provide the fi function for computing the 
interaction measure, and to provide the mapping from the lattice variables to 6/(/^ and aniq 
to the physical variables T and m^/mp. At the lattice spacing used here, lattice effects on 
the thermodynamics are large and flavor symmetry is strongly broken. We expect to pursue 
this project at smaller lattice spacing in the future. 
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TABLE I. Spectrum runs used to construct the nonperturbative QCD (3 function. We also 
tabulate the vr and p masses, together with the values of our interpolating expression at these 
points, using the functions in Table III. Note that the points with = 5.7 were not used in the 
fitting. 



run 


ref. 


aniq 




lattice size 


m„ 


nip 






1 


[9] 


0.1 


5.35 


8^ X 24^ 


0.8019(1) 


1.432(5) 


0.805 


1.438 


2 


[10] 


0.1 


5.375 


8^ X 24 


0.8088(7) 


1.408(11) 


0.811 


1.399 


3 


[11] 


0.1 


5.4 


6^ X 24 


0.819(4) 


1.38(1) 


0.814 


1.362 


4 


[10] 


0.1 


5.525 


8^ X 24 


0.8262(8) 


1.210(16) 


0.817 


1.211 


5 


[10] 


0.05 


5.32 


8^ X 24 


0.5827(7) 


1.406(48) 


0.543 


1.315 


6 


[15] 


0.05 


5.47 


8^ X 24 


0.6112(5) 


1.023(2) 


0.611 


1.024 


7 


[10] 


0.025 


5.2875 


8^ X 24 


0.4184(9) 


1.342(46) 


0.375 


1.274 


8 


[12] 


0.025 


5.445 


16^ X 24 


0.4488(4) 


0.918(4) 


0.447 


0.912 


9 


[13] 


0.025 


5.6 


16^ X 32 


0.4197(7) 


0.6396(28) 


0.425 


0.641 


10 


r ^ a1 

[14] 


0.025 


5.7 


32-^ X 32 


0.383(1) 


0.530(2) 


0.390 


0.511 


11 


[16] 


0.02 


5.5 


16^ 


0.3901(17) 


0.758(48) 


0.403 


0.769 


12 


[16] 


0.02 


5.6 


16^ 


0.3696(31) 


0.563(11) 


0.377 


0.599 


13 


[16] 


0.02 


5.7 


20* 


0.3402(17) 


0.4916(30) 


0.338 


0.464 


14 


[12] 


0.0125 


5.415 


16^ X 24 


0.3239(5) 


0.883(6) 


0.330 


0.891 


15 


[16] 


0.01 


5.5 


164 


0.2942(44) 


0.636(37) 


0.293 


0.692 


16 


[13] 


0.01 


5.6 


16^ X 32 


0.2667(8) 


0.5133(22) 


0.263 


0.512 


17 


[16] 


0.01 


5.7 


20^ X 204 


0.2451(23) 


0.4184(70) 


0.219 


0.368 


18 


[17] 


0.004 


5.48 


16^ X 32 


0.189(1) 


0.62(4) 


0.197 


0.686 
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TABLE II. Nonperturbative QCD (3 function. The upper section is for the (3 function calculated 
by the direct method. The first two entries are from new simulations and correspond to the 
VT p and the PV p, respectively. In these new simulations the couplings were actually varied 
asymmetrically in space and time around the cited values. The other points are from the literature 
and correspond to the VT p. "runs" gives the runs from Table I used to construct the (3 function, 
"corr" is the scaled correlation between the two components of the (3 function. The lower section 
of the table contains the (3 function obtained from the fits to the mass spectrum in Table III at a 
few selected points. Finally, we quote the perturbative value. 





aniq 


runs 




daniq 
(91n(ot) 


dlii(amq) 
(91n(ot) 


corr. 


5.35 


0.1 


1 


-1.201(342) 


0.304(20) 


2.04(20) 


-0.99 


5.35 


0.1 


1 


-0.379(128) 


0.253(8) 


1.53(8) 


-0.97 


5.379 


0.075 


12 5 6 


-0.333(11) 


0.1774(9) 


1.37(1) 


-0.50 


5.381 


0.0375 


5 6 7 8 


-0.243(15) 


0.0886(8) 


1.36(2) 


-0.77 


5.3825 


0.0208 


7 8 14 


-0.251(23) 


0.0444(5) 


1.13(2) 


-0.84 


5.465 


0.0169 


8 11 14 15 


-0.176(53) 


0.0351(15) 


1.08(9) 


0.40 


5.505 


0.0333 


6 8 9 


-0.338(58) 


0.067(2) 


1.01(6) 


0.98 


5.55 


0.01625 


9 11 15 16 


-0.249(88) 


0.027(3) 


0.64(18) 


0.98 


5.65 


0.0175 


9 10 16 17 


-0.344(35) 


0.024(2) 


0.37(11) 


0.84 


5.35 


0.1 




-0.565(25) 


0.200(4) 


1.00(4) 


0.52 


5.35 


0.025 




-0.278(14) 


0.082(2) 


2.28(8) 


0.09 


5.55 


0.025 




-0.226(3) 


0.0450(2) 


0.80(8) 


0.49 


5.35 


0.0 




-0.371(20) 








5.55 


0.0 




-0.270(4) 








00 


0.0 




-0.734 
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TABLE III. Fits to spectrum data used to construct the nonperturbative QCD (3 function. 
(Not recommended for > 5.6) 

m^/mp = ^0777^(5. 10 + 12.89(6/^2 _ 5 45^ _ 2.49(6/^2 _ 5.45)2^ 
-0771,(15.05 + 34.38(6/^2 _ 5 45^^ _^ 16.51(amg f/'^ 
with 8 degrees of freedom 23.4. 
Confidence level 0.0029. 

Runs from Table I used in fit 1, 2, 3, 4, 5, 6, 7, 8, 9, 11, 14, 15, 16, 18 
nip = 0.72 - 2.25(6/^2 _ 5 45^ _^ 1.75(6/^^ _ 5 45^2 _^ 7.75^^^ 

+ 10.010777,(6/^2 _ 5 45^ _ 20.08(am,)2 
^2 with 8 degrees of freedom 20.7. 
Confidence level 0.008. 

Runs from Table I used in fit 1, 2, 3, 4, 5, 6, 7, 8, 9, 11, 14, 15, 16, 18 
= 0-49 - 2.20(6/^2 _ 5 45^ _^ 3.00(6/^2 _ 5 45^2 _^ ii.02am, 
+ 7.320777,(6/^2 _ 5 45^ _ 36.61(am,)2 
^2 with 5 degrees of freedom 8.9. 
Confidence level 0.11. 

Runs from Table I used in fit 1,4,8,9,11,12,13,14,15,16,18 
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5.2 5.3 5.4 5.5 5.6 

6/g^ 

FIG. 8. The effect of tfie finite step size on tlie pressure. Tlie fancy plusses are from step size 
0.03, tlie fancy crosses from 0.02. Tlie octagons give the result extrapolated to zero step size. The 
bursts show the mass integration results. 
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6 




5.2 5.4 5.6 5.8 

6/g^ 

FIG. 9. The pressure in units of T"^ as a function of gauge coupling The bursts are 

extrapolations of the mass integrations to anig = (see Fig. 10 and Eq. 26). 
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6 




1 2 



rriq/T 

FIG. 10. The pressure in units of T"* as afunction of the quark mass nig/T from the integrations 
over mass. The diamonds show give = 5.53, the octagons = 5.45 and the squares the 

= 5.35 results. The bursts are extrapolations of the mass integrations to anig = (see Eq. 26). 
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5.2 5.4 5.6 5.8 




FIG. 11. The energy density in units of T"^ as a function of the gauge coupling. The lower 
curves are the pressure values displayed for comparison. The octagons are for anig = 0.025, or 
niq = T/10, and the squares are for ania = 0.1, or niq = OAT. The errors contain the uncertainty 
in the (3 function. The anig = 0.1 energy could not be computed at higher because of the lack 
of a reliable /3 function in this region of and anig. 
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0.0 0.1 0.2 0.3 0.4 

T(GeV) 

FIG. 13. The energy density in units of T"^ as a function of the temperature. The lower curves 
are three times the pressure displayed for comparison. The octagons are for anig = 0.025, or 
niq = T/10, and the squares are for ania = 0.1, or nig = T/4. The errors contain the uncertainty 
in the /3 function. The arriq = 0.1 energy could not be computed at higher because of the lack 
of a reliable (3 function in this region of and anig. (The pressure curve for this case lies almost 
on top of the energy curve for anig = 0.025 and is easily overlooked.) The bursts are extrapolations 
of the mass integrations to nig = 0. 



32 



20 



15 



A 10 

CO 

I 



~i — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — r 



^ amq=0.1 



[] 
o 



$ $ ^ 



t] 





I ^ ^ I I I I I \ I I I I I I \ I I u 





5.30 5.35 5.40 5.45 



20 



15 



10 



~i — I — I — I — I — I — I — I — I — \ — I — I — I — I — I — I — I — I — r 



am =0.025 



cP 



]^-^ CD 

_l I \ I I I I I I I I I iCtl I L 



5.2 5.3 5.4 5.5 

6/g^ 



5.6 



FIG. 14. Gauge and fermion contributions to the interaction measure for (a) m^a = 0.1 and 
(b) niqa = 0.025. The circles give the gauge part, and the diamonds give the fermion part. The 
squares are the total. 
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